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ABSTRACT 

An explanation of long-lived Saturn’s North Polar hexagonal circumpolar jet 
in terms of instability of the coupled system polar vortex - circumpolar jet is 
proposed in the framework of the rotating shallow water model, where 
scarcely known vertical structure of the Saturn’s atmosphere is averaged out. 
The absence of a hexagonal structure at Saturn’s South Pole is explained 
similarly. By using the latest state-of-the-art observed winds in Saturn’s 
polar regions a detailed linear stability analysis of the circumpolar jet is 
performed (i) excluding (“jet-only” configuration), and (2) including 
(“jet+vortex” configuration) the north polar vortex in the system. A domain 
of parameters: latitude of the circumpolar jet and curvature of its azimuthal 
velocity profile, where the most unstable mode of the system has azimuthal 
wavenumber 6, is identified. Fully nonlinear simulations are then performed, 
initialized either with the most unstable mode of small amplitude, or with 
the random combination of unstable modes. It is shown that developing 
barotropic instability of the “jet+vortex” system produces a long-living 
structure akin to the observed hexagon, which is not the case of the 
“jet-only” system, which was studied in this context in a number of papers 
in literature. The north polar vortex, thus, plays a decisive dynamical role. 
The influence of moist convection, which was recently suggested to be at the 
origin of Saturn’s north polar vortex system in the literature, is investigated 
in the framework of the model and does not alter the conclusions. 
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1 Introduction 


Visible imagery acquired during the Voyager 1 & 2 fly-by of Saturn in 1980-1981 
revealed the hexagonal cloud pattern of Saturn’s north pole at about 77°N plane- 
tographic latitude (Godfrey, 1988). As followed from these early observations, the 
motion of the associated cloud structures suggested that the structure is related to 
a circumpolar jet-stream, but the hexagon pattern itself appeared to be stationary 
(i.e. very close to the rotation of Saturn’s interior). The early diagnostics were 
confirmed first in the 1990s by Hubble Space Telescope observations (Sanchez- 
Lavega et al., 1993), and then in the 2000s by the Cassini orbiter firstly in thermal 
infrared images (Baines et al., 2009), followed by visible imagery once Saturn’s 
northern pole came out of the winter darkness. The high spatial resolution of the 
Cassini visible and infrared images permitted to evidence the presence of the north 
polar vortex (NPV), in addition to the hexagonal jet. The exceptional duration 
of the Cassini mission allowed Sanchez-Lavega et al. (2014) to conclude that the 
hexagon resists seasonal changes, before Antunano et al. (2015) compiled the re- 
peated cloud-layer observations of the hexagon to obtain a complete profiling of 
the winds in this structure. Thus far, such hexagonal feature has not been observed 
at Saturn’s South pole. 

The persisting hexagonal pattern in polar projected maps demonstrates that a 
prominent wavenumber 6 perturbation shapes Saturn’s circumpolar jet stream at 
latitude ~77°N. An early interpretation by Allison et al. (1990) was that the an- 
ticyclonic North Polar Spot (NPS) vortex (not to confuse with the polar vortex), 
visible at that time in the vicinity of the hexagon forced a stationary planetary 
(Rossby) wave which gave the polar jet its hexagonal shape, but the apparent 
disappearance of the NPS between 1995 and 2004 invalidated this scenario (see 
Sayanagi et al. (2016) for the history of the observations of the NPS and the 
hexagon). The nature and origin of the hexagon has been addressed since then 
both in laboratory experiments and numerical models. Polygonal structures are 
often reported in rotating flow experiments, although the wavenumber 6 config- 
uration remained elusive (Vatistas et al., 1994; Marcus and Lee, 1998; Jansson 
et al., 2006; Bergmann et al., 2011). Aguilar et al. (2010) discussed more specif- 
ically the relevance of rotating tank experiments to modeling Saturn’s polar jet 
and concluded that the hexagonal structure could be resulting from the barotropic 
instability of the jet, and that the wavenumber 6 prominence is conditioned by the 
intensity of the jet and bottom friction. They were first to notice a dependence 
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of the azimuthal wavenumber of the most unstable mode on the deformation ra- 
dius in their instability analysis in the framework of the quasi-geostrophic (QG) 
model. Using a general circulation model with a domain limited to Saturn’s north- 
ern polar region, and a polar jet as initial condition, Morales-Juberias et al. (2011) 
reproduced the structures found by Aguilar et al. (2010) in laboratory experi- 
ments: the hexagonal shape resulting from a vortex street formed by developing 
barotropic instability of the jet, with cyclonic and anticy clonic vortices at poleward 
and equator- ward sides of the jet, respectively. Yet, the vortices forming the street 
were too large and too strong, with respect to observations, and this mechanism 
was giving rise to larger propagation speed than the observed one (Sanchez-Lavega 
et al., 2014). Morales-Juberias et al. (2015) proposed an alternative “meandering 
jet” model which matches the morphology (including sharp potential vorticity, 
hereafter PV, fronts) and phase speed of Saturn’s hexagon, provided that an ad 
hoc vertical shear of the jet is introduced. The same authors concluded that deep 
jets evolve into vortex streets and shallow jets evolve into meanders, both being 
possibly at the origin of the absence of seasonal variability of the hexagonal jet. 
Although the meandering jet model of Morales-Juberias et al. (2015) offers thus 
far the best agreement with the characteristics of the observed Saturn’s hexagon, 
the putative meridional temperature gradients accompanying the vertical shear of 
the (supposedly shallow) polar jet are yet to be confirmed. 

Despite these developments, there is not enough convincing dynamical explana- 
tion of the existence and origin of Saturn’s North polar hexagon, as well as the 
absence of its counterpart at Saturn’s South pole. In the present paper we pro- 
pose an alternative approach to the problem. Our analysis is performed in the 
framework of a simple rotating shallow water (RSW) model, which is of use for 
modeling atmospheres of giant planets, cf. (Dowling, 1995). A specificity of our 
approach is that the model is considered in the polar tangent plane approxima- 
tion, the so-called gamma-plane, in order to take into account the variations of the 
Coriolis parameter with latitude in the polar regions. The use of the model, which 
may be obtained by vertical averaging of the full primitive equations, is justified 
by the fact that information on the vertical structure of Saturn’s atmosphere is 
scarce. All vertical structure being averaged out, the model is barotropic, in the 
sense that it does not contain any vertical shear, although it does allow for vertical 
stretching of fluid columns. In the same sense the instabilities in the model will 
be called barotropic below. The model combines simplicity with dynamical con- 
sistency and, with the astronomical parameters of Saturn being given, contains a 
single adjustable parameter: the effective Rossby deformation radius, Rd ■ The QG 
model used for theoretical analysis in Aguilar et al. (2010) can be obtained from 
RSW in the limit of vanishing Rossby numbers. The RSW model, which allows for 
efficient high-resolution numerical implementations, has been used for the analysis 
of various features of the Saturn’s atmosphere: mid-latitude jets (Showman, 2007), 
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equatorial super-rotating jet (Scott and Polvani, 2008), NPV (O’Neill et al., 2016) 
(in two-layer version with a moist-convective forcing in the latter work). Yet, to our 
knowledge, there are no studies in the literature applying the RSW model to the 
developing barotropic instability of the circumpolar Saturn’s jet. The particularity 
of our approach is that we explore the dynamics in two different configurations 
of the mean zonal velocity profile: (i) considering only the 77°N circumpolar jet, 
with observed velocity profile, and no central vortex at 90°N (i.e. the traditional 
“jet-only” configuration considered in all above-cited papers), and (ii) considering 
the full “jet + vortex” system with the observed zonal velocity profile, which was 
not considered in the existing literature in this context. The state-of-the-art mean 
zonal velocity profiles on Saturn (Antunano et ah, 2015) are used in linear stability 
studies and for initializations of nonlinear simulations. 

The paper is organized as follows: in section 2 we introduce the model, discuss its 
general features and vortex solutions, present the analytical fits of the observed 
velocity profiles, and formulate the linear stability problem. Section 3 is devoted to 
the study of the jet-only configuration, with presentation of the results of the linear 
stability analysis and nonlinear simulations of the saturation of the instability. In 
section 4 we analyse the jet+vortex configuration along the same lines. Finally, in 
section 5 we present our conclusions and discussion. A discussion of the influence 
of the effects of moist convection upon the evolution of the instability is given in 
the Appendix. 


2 Rotating shallow water model for large-scale atmospheric motions, 
vortex solutions and their (in)stability 


2.1 The model 


Rotating shallow water equations, in the absence of forcing and dissipation read: 

7^7 + (/ + ~)(z x v) = -gVh, (2.1) 

1)1 r 

Dh 

+ = 0. (2.2) 

Here, by anticipating the application to polar jets and vortices, we use the cylin- 
drical coordinates on the plane (r, 9), and z is the unit vector normal to the plane. 
v = ur + v6 is the horizontal velocity with u. v denoting radial and azimuthal 
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components, respectively, D / Dt = d/dt + v • V denotes the material derivative, 
V = fd r + 9(l/r)de is the gradient in the plane, where d r and 3q denote partial 
derivatives with respect to corresponding arguments, h is the thickness of the layer, 
g is gravity and / is the Coriolis parameter. 

Bottom topography b(r, 9) can be easily introduced in the model by replacing h by 
h — b in (2.2), and can be used to parameterize the deep layers of the atmosphere. 
This would however introduce ad hoc parameters, which we want to avoid. 

The RSW equations possess a Lagrangian invariant, the potential vorticity (PV): 

z • (V x v) + / 

9 = h ’ 

where z ■ (V x v) is relative vorticity. The PV anomaly 

PV A = q — /o/ Hq, (2.4) 

where f 0 and H 0 are reference values of Coriolis parameter and thickness of the 
layer (see below), will be used for diagnostics in what follows. We consider the 
RSW equations in the polar tangent plane, the so-called y-plane, where the effects 
of sphericity in the Coriolis parameter are taken into account in the lowest-order 
approximation: 

f = fo sin(<#>) = /u (l - + ■■■) = fo (l - 7‘f* 2 + •••) ; 0 < r* < A. 

(2.5) 

Here f 0 = 2 fl s , / is planetographic latitude, and R s represent Saturn’s angular 
velocity and radius, respectively (we will neglect the effects of non-sphericity in 
what follows). Approximate values for the mean Saturn radius and fo are respec- 
tively 55000 km and 3.2 x 10 -4 s _1 (Read et al., 2009b; Sanchez-Lavega et ah, 
2014). As is well-known, the rotating shallow water system possesses an internal 
scale, the Rossby deformation radius R r [ = sfgHo / fo . 

Equations (2.1), (2.2) can be obtained by vertical averaging between two material 
surfaces of the full three-dimensional primitive equations for the atmosphere in 
pseudo- height isobaric coordinates, e.g. (Bouchut et al., 2009), when one of these 
surfaces is considered fixed (a constant pressure level), and one is free. The poten- 
tial temperature is considered uniform through the layer in this approximation, 
although an extension of the model including horizontal potential temperature 
gradients is possible. In this case g is the gravity acceleration of Saturn, and H 0 is 
the thickness of the layer. The model can be also considered as a two-layer model 
with strong disparity in depths and/or densities of the layers. In this case H 0 plays 
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a role of equivalent depth, and g is so-called reduced gravity, i.e. gravity weighted 
with the stratification parameter. The deformation radius in this interpretation 
is related to stratification properties of the atmosphere. Whatever the interpreta- 
tion is, at a given rotation rate /o/2 , the deformation radius R,j is the only free 
parameter of the model, as the velocity scale can be taken to be \/gH () . 

In terms of velocity components, the momentum and mass conservation equations 
of the model read: 


j^r ~ ~ ~ jv = - gd r h , 

Dt r 

(2,6) 

Dv uv 

(2.7) 

7T7 + ~+fu = -gdeh, 

Dt r 

d t h -\ — d r (hru ) -I — de(hv) = 0. 

(2,8) 


2.2 Vortex solutions 


Stationary axisymmetric vortex solutions of the equations (2.6) - (2.8) are those 
with zero radial velocity u = 0, some axisymmetric distribution of azimuthal ve- 
locity v = V(r), and the thickness field H(r ) related to V(r) through the cyclo- 
geostrophic (gradient wind) balance: 

— + fV = gd r H. (2.9) 

r 

For a given V(r), corresponding thickness profile H(r) can be obtained by inte- 
gration of equation (2.9). Note that in the QG model used for theoretical analysis 
in Aguilar et al. (2010) the centrifugal acceleration (the first term in (2.9)) is ne- 
glected. While this approximation is fully justified for the circumpolar jet, it is 
less so for the central vortex. Herein we consider solutions of (2.9) with V(r ) cor- 
responding to the observed time- and zonally averaged profiles of Saturn’s zonal 
wind (Antunano et ah, 2015). These profiles, with error margins, are represented 
by dashed lines in Figure 1. 

It is useful to have a simple analytic expression mimicking the observed velocity 
profiles. Although the linear stability analysis can be accomplished directly with 
digitized observational velocity profile, it is convenient to control the shape of 
the velocity profile in terms of a small number of parameters. Dependence on 
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Figure 1. Observed zonal velocity profiles in the northern (left panel) and southern 
(right panel) hemispheres and their fits used in the paper. Three dashed lines represent 
the observed non-dimensional mean and error margins of the averaged zonal velocity 
(Antunano et al., 2015), and the solid line is the composite of two model profiles (in 
red) fitting the NPV and the circumpolar jet. The fit of the velocity profile of the 
NPV gives e — 0.15, a — 0.42, fd — 1.3, ro = 0,m = 1 and that of the circumpolar jet 
gives e = 0.08, a — 0,/3 = 2,ro = 3.37, m — 3. Similarly, south polar vortex corre- 
spond to e = 0.16, a = 0.5 , (3 = 1.2, ro = 0,m = 1 and southern circumpolar jet jet to 
e — 0.08, a — f), fd — 2,r^ — 4.6, rn = 3. R ( i ^ 3200 km is taken for both poles. 




Figure 2. PV distributions corresponding to observed velocity profiles of Fig. 1 in the 
northern (left panel) and southern (right panel) hemispheres and their analytic fits. 

these parameters could be then analysed. The PV, which is a key characteristic of 
the flow, is very noisy, if determined directly from the observed profiles, cf. Fig. 
2, whereas a PV field computed from the analytic fit enables a straightforward 
analysis of possible instability, see below. For initialization of the non-linear RSW 
simulations, the use of the analytic fit allows for straightforward integration of 
(2.9), and thus diminish the discretisation errors. Therefore, we fit the peaks of 
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the velocity distribution with the help of the following simple formula: 

V(r) = e(r — ro) a e~ rn ^ r ~ r °' >13 , a,/3,e,m>0, (2-10) 

where e measures the intensity of the velocity held, ro tunes the distance of the 
velocity peak of the jet from the pole, and other parameters allow to fit the shape 
of the distribution (r 0 = 0 for the central vortex and a = 0 for the jet). All velocity 
profiles are normalized by their maximum values, so the maximum value of velocity 
is equal to e. By applying this formula to both central vortex and circumpolar jet 
(red curves in Fig. 1), we get synthetic profiles represented by solid black lines in 
the figure. We do not seek to reproduce the “shoulder” in the velocity profile of 
the central vortex, as it turns out that it is inessential for the dominant long-wave 
instability of the system which we are interested in (a test study was performed 
to check this fact). 

The stability of any velocity profile V (r) (and corresponding H (r) obtained from 
the cyclo-geostrophic balance) with respect to small perturbations can be analyzed 
by standard means. We represent each field as the solution in question and a small 
perturbation (denoted by a prime): 

u(r,9,t) — u'(r,9,t), (2-11) 

v(r, 9, t) — V(r) + v'(r, 9, t), (2-12) 

h(r,9,t) — H(r) + h!(r,9,t), (2-13) 

and inject these expressions in the full system (2.6) - (2.8). Retaining only the linear 
contributions in perturbation terms, we are looking for normal-mode solutions with 
harmonic dependence on time and polar angle 

(■ u' ,v' ,h!)(r,9,t ) = Re[(iu, v, h){r)e l ^ ie ~ ut ^]^ (2-14) 

where / and co are the azimuthal wavenumber and the complex eigenfrequency, 
co = Ur + iuJi. We thus formulate the problem as an eigenvalue problem for eigen- 
frequencies co. A positive imaginary part for an eigenfrequency coj corresponds to 
an instability with a linear growth rate a = coj. 

Following common practice, it is convenient to formulate this eigenproblem in non- 
dimensional terms. By using the deformation radius Rj as the scale for r, f ( J 1 as 
the scale for t, and \fgHo as the scale for velocity, the eigenproblem in question 
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(2.15) 


where variables marked with the asterisks are non-dimensional, and D r * denotes 
the differentiation with respect to r*. The non-dimensional 7* = O(10 -3 ). It should 
be noted that the parameter e acquires in this way a meaning of the Rossby number. 


Pseudospectral collocation method (Trefethen, 2000) is used for linear stability 
analysis of the eigen- problem (2.15). The system is discretized over A-point grid 
and D r * becomes the Chebyshev differentiation operator. To avoid the Runge 
phenomenon, Chebyshev collocation points are used, following Lahaye and Zeitlin 
(2015) where the method of Boyd (1987) was adapted to the stability problem of 
circular vortices. 


3 Instability of the “jet-only” configuration and its nonlinear satura- 
tion 


In this section, we analyse the instability of the North Polar circumpolar jet in 
the absence of the central vortex (“jet-only” configuration), following the existing 
studies in the literature (Aguilar et ah, 2010; Morales- Juberias et al., 2011, 2015), 
albeit in a different model. The velocity profile of the jet is fitted by a Gaussian, 
i.e. the profile (2.10) with a = 0. The non-dimensional amplitude of the jet is 
e = 0.08. The values of parameters m and r 0 corresponding to the most reasonable 
fit of the observations with error bars (see Fig. 1) are m = 3 and r 0 = 3.37. 


3. 1 Results of the linear stability analysis 


We present in this subsection the results of linear stability analysis for the “jet- 
only” configuration, carried out along the lines of section 2. The barotropic in- 
stability of zonal jets, as well as its nonlinear saturation, are well understood in 
geophysical fluid dynamics. In the framework of the RSW model sufficient condi- 
tions of stability of zonal jets are given by Ripa’s criteria (Ripa, 1990). Like for the 
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famous Rayleigh - Kuo stability criterion (which should be adapted, however, to 
the current 7- plane context) the meaning of this criterion is that in order to have 
the instability the PV gradients should change sign within the flow. Physically, 
this means that there should be a possibility for Rossby waves, which owe their 
existence to PV gradients, to propagate in opposite directions in the flow, phase- 
lock, and grow in amplitude. As follows from Fig. 2 there is indeed a change of 
sign of the gradient of PV, which is by the way more pronounced in the Northern 
hemisphere, and hence we expect the barotropic instability. The main question 
then is: under what precise circumstances the mode with azimuthal wavenumber 
l = 6 is the dominant instability mode of the circumpolar jet? Our linear stability 
analysis demonstrates that for the location of the jet at the distance from the 
pole r 0 corresponding to the observations, the dominant instability mode has an 
azimuthal wavenumber 1 = 6 for a choice of Rj ~ 3200 km. This value of Rd is in 
line with the values of 2500 — 3000 km adopted in Aguilar et al. (2010) and O’Neill 
et al. (2016) (cf. their appendix 2), and are of the same order of magnitude as the 
estimates from observations in Read et al. (2009a), although about twice larger. 
Yet Rd remains to be better constrained by observations. 

Changing the radial (latitudinal) position of the jet, 77, changes the wavenumber of 
the most unstable mode. This is shown in Fig. 3, where the growth rates of unstable 
modes as a function of e and tq are displayed in a range of values of e and tq covered 
by our linear stability analysis). Fig. 3 shows that the azimuthal wavenumber of 
the most unstable mode varies from l = 2 to l = 8 with r 0 increasing from 2 to 5 
in non-dimensional terms. For some values of r 0 , several modes exhibit very close 
growth rates, especially at small e: in such configuration a combination of these 
modes would eventually shape a circumpolar jet without pronounced polygonal 
pattern. It is worth noting that an increase of velocity amplitude of the jet, e, does 
not change the azimuthal wavenumber of the most unstable mode, but increases 
the growth rate, cf. Fig. 3. 

It should be emphasized, however, that even a small difference in growth rates 
yields exponentially- growing differences between corresponding wave patterns, whence 
the primary importance of the leading mode. 

In Figure 4 we present the regions in the parameter plane r$ — m where unstable 
modes with a given l have the highest growth rate at fixed Rd . Higher values of m 
correspond to stronger curvatures of the velocity distribution. The value of Rd also 
influences the results. For larger values of Rd, the observed position of the maxi- 
mum jet velocity shifts closer to the pole in non-dimensional terms and, as follows 
from the results described in the previous paragraph, the azimuthal wavenumber of 
the leading instability becomes smaller, and vice-versa. By reversing the argument, 
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Figure 3. Dependence of the growth rate of the unstable modes on azimuthal 
wavenumber for different latitudinal positions of the jet. V(r) is Gaussian , with 
m = 3, ol — 0, (3 — 2, Rd « 3200 km. 


from our analysis and from the fact that the observed dominant wavenumber is 
l — 6, we can infer an estimate for R & « 3200 km. 


An important unsolved problem is the difference in morphology between the hexag- 
onal northern polar jet and the southern polar jet which does not exhibit a pro- 
nounced polygonal form in Saturn’s atmosphere, as follows from the observations 
(Antunano et ah, 2015). Our linear stability analysis of “jet-only” configuration 
provides a clue for explanation. As follows from Fig. 1, at a given value of the 
value of r 0 is higher in the southern hemisphere (i.e. the southern jet is farther 
from the pole), which yields a leading instability with l > 8 (Fig. 3), and hence no 
hexagonal pattern at Saturn’s South Pole, although a partial polygonal structure 
was discussed in the literature (Vasavada et ah, 2006). 
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Figure 4. Left panel: Distribution of azimuthal wavenumbers of the most unstable modes 
in the plane of parameters vq — m. Black (white) square corresponds to the observed pro- 
file in the northern ( southern ) hemisphere. V (r) is Gaussian, with e = 0.08, a = 0, (3 = 2. 
Azimuthal wave-numbers larger than eight are not included. Right panel: Dependence 
of the growth rates of the unstable modes of the “jet-only” configuration in northern 
hemisphere of the value of equivalent depth. Variation of the growth rate with azimuthal 
wavenumber for the jet corresponding to the black square in the left panel corresponds 
to Rd ~ 3200 km. 


3.2 Non-linear saturation of the instability 


In order to study the nonlinear saturation of the instability identified in the pre- 
vious subsection, we perform direct numerical simulations with a well-balanced 
entropy-satisfying finite-volume scheme (Bouchut, 2007) for the RSW equations. 
The numerical scheme was successfully tested in studies of the saturation of the 
barotropic instability of jets and vortices (Lambaerts et ah, 2011; Lahaye and 
Zeitlin, 2015). The scheme has no explicit dissipation, only a numerical one. The 
dissipation is concentrated in the zones of strong gradients of velocity and thick- 
ness. For flows with low Rossby numbers, as is the case here, the scheme is prac- 
tically dissipationless and allows for high-resolution long-time simulations. The 
simulations below are initialized with the most unstable mode of small amplitude 
(several per cent with respect to the background flow) superimposed onto the jet, 
or with an ensemble of unstable modes with random phases. 

The evolution of the instability follows the well-known scenario of the saturation 
of the barotropic instability. Namely, a series of vortices of opposite sign appear 
at the crests and lows of the initial unstable wave, cf. Fig. 5. This is similar to the 
vortex-street patterns observed in simulations of Morales-Juberias et al. (2011) 
and laboratory experiments of Aguilar et al. (2010). As time goes on, the vor- 
tices intensify, resulting in a strong deformation of the polar jet. Although the 
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Figure 5. Snapshots of the pressure (left panel) and potential vorticity anomaly (right 
panel) at t = 150/o -1 during the saturation of the instability in the “jet-only” configu- 
ration. 


outcome does preserve a six-fold symmetry of the initial unstable mode, the ob- 
served hexagon with quasi-straight sides (corresponding to sharp PV gradients) 
is not reproduced. Moreover, as follows from Fig. 6, the azimuthally averaged ve- 
locity of the resulting structure displays formation of counter-jets at both sides 
of the original jet, which is typical for the saturation of the barotropic instabil- 
ity (e.g. (Lambaerts et al., 2011)), and does not match the observations. At later 
stages the hexagonal structure practically disappears (cf Fig. 7). The evolution 
of radial distribution of zonally averaged PVA, as follows from the simulations, is 
presented in Fig. 6. As follows from the figure, the flow reorganizes by diminishing 
PV gradients and flattening the radial PV distribution in the region between the 
circumpolar jet and the NPV , and thus tending to “cure” the instability, but this 
process is not yet over at t — 450/q” 1 . Indeed, opposite PV gradients of comparable 
strength are still present in the flow at this stage, thus maintaining the barotropic 
instability. We do not trace further the process of saturation which is well-known 
and leads to complete “healing” of the instability, appearance of counter-jets (al- 
though here we have a new element with respect to classical studies element, the 
7- effect), and does not produce a quasi-stationary hexagonal structure. Note that 
the (diminishing) maximum of PVA shifts poleward during the saturation of the 
instability. 

Although the l = 6 instability is dominant, the growth rates of l = 5 and l = 7 
modes are rather close to the growth rate of l = 6 (Fig. 3, left bottom panel). If 
the initialization of direct numerical simulations is made with a combination of 
these modes with equal amplitudes and random phases, the velocity distribution 
is still meandering and possesses zonal counter-jets at the intermediate stages, and 
does not exhibit the hexagonal shape at late stages (not shown). 
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Figure 6. Left panel: Superposition of the observed mean azimuthal velocity profile with 
error margins (thin dashed), model azimuthal velocity profile (solid), and mean az- 
imuthal velocity profile (dashed) at t = 150/o -1 during the saturation of the instability 
in the “jet-only” configuration. Right panel: Evolution of the radial distribution of zon- 
ally averaged PVA in the ‘ jet-only ” configuration. 



Figure 7. Same as in the left panel of Fig 5, but for 
t = 350/J" 1 , (left panel), and 450 f^ 1 (right panel). 


4 Instability of the “jet + vortex” configuration and its nonlinear sat- 
uration 


It should be stressed that the “jet-only” approximation of section 3, which is 
being overwhelmingly used in the literature thus far, would be valid if the jet 
and the vortex were spatially well separated. As follows from Fig. 1, the velocity 
peaks are separated by approximately 3 Rd- A well-known property of the RSW 
is that the deformation radius Rd plays a role of a “screening” radius beyond 
which the velocity field generated by a vortex falls off exponentially. In spite of 
their rather well-separated velocity peaks, the velocity distributions of the central 
vortex and of the circumpolar jet are not sufficiently far from each other to prevent 
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Re(ro)= 0.044842, lm(co)= 0.01 7369 


Re(co)= 0.049225, lm(©)= 0.01 6855 




Figure 8. Structure of the most unstable mode with l = 6 of “jet-only” configuration (left 
panel) and of “jet-h vortex” one (right panel)in the northern hemisphere corresponding 
to the background velocity profile of Fig.l. Vertical lines represent positions of critical 
levels. 

the central vortex influencing the perturbations of the jet. We will show below that 
this influence changes the evolution of the instability, although the most unstable 
mode of the composite jet + vortex system remains practically the same as in the 
“jet-only” case. 


4-1 Results of the linear stability analysis 


The most important result following from the linear stability analysis in the 
“jet+ vortex” configuration (Fig. 8) is that the most unstable mode has l = 6 
as in the “jet-only” case, and that its structure is practically the same. The big 
picture of the dependence of the azimuthal wavenumber of the most unstable mode 
on the parameter ro is also close to what was observed in the “jet-only” configu- 
ration, as follows from Fig. 10. The configuration corresponding to the South Pole 
has the most unstable mode with l = 8 (at the value R 4 ~ 3200 km used as a 
reference for comparison with the “jet-only” case). The position of the southern 
polar jet being much farther from the pole in comparison with the northern polar 
jet, we can expect, as was the case in the “jet-only” configuration, that instability 
happens at higher l than the l = 6 leading to the hexagonal structure. The most 
unstable l can be pushed even higher by diminishing R,i, as follows from Fig. 11. 
The configuration corresponding to the South Pole has the most unstable mode 
with l = 8 (at the value R r i ~ 3200 km used as a reference for comparison with the 
“jet-only” case). The position of the southern polar jet being much farther from 
the pole in comparison with the northern polar jet, we can expect, as was the case 
in the “jet-only” configuration, that instability happens at higher l than the l = 6 
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Figure 9. Same as in Fig. 4, but for the “jet-hvortex” configuration. The white square 
corresponds to the parameters in the right panel of Fig. 1 . 
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Figure 10. Same as in Fig. 4, but for the “jet+vortex” configuration. The white square 
corresponds to the parameters in the right panel of Fig. 1 . 

leading to the hexagonal structure. The most unstable l can be pushed even higher 
by diminishing R ( [, as follows from Fig. 11. As one can infer from the observations 
(Antunano et ah, 2015) the southern circumpolar jet has a much higher azimuthal 
wavenumber than l = 8 (l = 10 — 11), we thus may conclude that the equivalent 
depth near the South Pole is, by some reason, smaller than in the vicinity of the 
North Pole. 


4-2 Non-linear saturation of the instability 


Direct numerical simulations of the evolution of the instability of the “jet+vortex” 
system were performed using the same procedure as in the “jet-only” case. We 
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Figure 11. Linear growth rates o of different azimuthal modes of “jet+vortex” configu- 
ration at the North Pole (left panel) and South Pole (right panel) at different values for 
R,]. The non-dimensional velocity profiles are the same as Fig.l. 

present in Fig. 12 the snapshots of pressure and potential vorticity anomaly at 
times t = 150/o -1 and t = lOOO/o -1 clearly showing the formation and mainte- 
nance of a hexagonal structure matching the observed feature. Compared with the 
evolution of the “jet-only” configuration in Figs. 5, 7 the “jet+vortex” case ex- 
hibits much less deformation of the jet with time, resulting in the long-living quasi- 
stationary hexagonal pattern exhibiting sharp straight PV boundaries, which was 
not the case in the “jet-only” configuration. The evolution of the zonally averaged 
PVA is presented in Fig. 13. As follows from the figure, the jet+vortex configura- 
tion “cures” the instability at t = 1000/T 1 . Positive PV gradients become small 
and are shifted farther from the zone of negative PV gradients of the jet. This pre- 
vents effective coupling of Rossby waves moving in opposite direction and, thus, 
eliminates the cause of instability. We performed a linear stability analysis of the 
mean zonal velocity profile corresponding to the PVA profile at t = 1000/q" 1 of 
Fig. 13 and found that although the unstable modes still exist, their growth rates 
are at least one order of magnitude less than those of unstable modes of initial 
profile. Hence, in accordance with qualitative considerations above, the instability 
does cure itself. This process, is, however, not sufficient to explain by itself why 
the hexagonal structure of finite amplitude persists. 

The inclusion of the NPV thus completely changes the evolution of the circumpolar 
jet. This could be further demonstrated by comparing the evolution of the mean 
azimuthal velocity in the “jet+vortex” case presented in Fig. 14 with that of the 
“jet-only” case in Fig. 6. The results are presented for two different intensities 
of the central vortex: one corresponding to the observations, and another one of 
the same shape but of higher amplitude. In the second case the system reaches a 
quasi-equilibrium configuration close to the observed one in high latitudes. Yet, 
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Figure 12. Upper row: Snapshots of pressure distribution in “jet+vortex” configuration 
at time = 150, 1000/o -1 respectively, left and right panels. Lower row: corresponding 
potential vorticity anomaly. V (r) is the same as Fig.l. For better visibility the colorbars 
are not the same in two lower figures. 

a discrepancy is observed in low latitudes, where the velocity field in simulations 
does not match observations. It should be stressed that our goal is, in the first 
place, to understand the dynamical role of the polar vortex, so we did not seek 
to represent the details of the velocity field beyond the jet with our simple fit. 
There are also technical reasons for this decision. First, the 7 plane approximation 
is losing its accuracy while moving farther from the pole, so trying to reproduce 
fine details of the velocity profile does not make much sense. Second, in order to 
prevent emitted inertia-gravity waves from returning back to the vortex, sponge 
boundary conditions were imposed at the boundary of the computational domain 
at r ~ 10. The mean velocity field, therefore should fall off to zero sufficiently 
far from the sponges, to be not perturbed by these latter. However, the (weak) 
retrograde jet which is present in the data could play a role. For example, an 
exchange of angular momentum between retrograde and prograde jets in Saturn’s 
atmosphere was reported in Liu and Schneider (2015). Simulations on the whole 
sphere are needed in order to correctly reproduce this part of the velocity profile, 
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Figure 13. Evolution of the radial distribution of zonally averaged PVA in the 
u jet+vortex ” configuration. 




Figure 14. Superposition of the observed zonal velocity profile with error margins (thin 
dashed), the model zonal velocity profile used in initialization of numerical simulations 
(solid), and zonal velocity profile after the hexagon formation (dashed). In the right 
panel the amplitude of the initial vortex was taken deliberately higher (e = 0.2) than in 
the right panel (e — 0.15), in order to arrive to a close to observations end state. 


which is beyond the scope of this paper. 

The same tendency is observed in the distribution of zonally averaged angular 
momentum represented in Fig. 15 by its deviation from the initial distribution. 
The angular momentum of the system is practically conserved (note that we have 
numerical dissipation and sponges at the boundaries). Fig. 15 gives an idea of 
angular momentum fluxes which are represented with the help of the difference of 
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Time = 1 00 to 200 [f” 1 ] Time = 1 00 to 450 [f” 1 ] 



Figure 15. Difference of zonally averaged angular momentum density per unit mass M of 
“jet” and “jet+vortex” configurations at times t — 200/q" 1 (left panel), and t — 450/q" 1 
(right panel), and time t — 100 f^ 1 , in the northern hemisphere. Vertical dashed line 
represents the initial maximum of of the jet velocity 


M(t) — M(£ 0 ) the angular momentum density per unit mass M 

r 2 

M(r)=rv + f — (4.1) 

at a given time t and a fixed time to- We have chosen t 0 = lfJO/, 7 1 to avoid 
oscillations due to emission of inertia-gravity waves, and their partial reflection by 
(imperfect) sponges, at the initial stages of the evolution. The Figure shows a net 
angular momentum flux from the polar vortex to the jet, which diminishes in time 
in the jet+vortex configuration, but continues to act in the jet-only one. 

A tentative interpretation of the stabilization of the hexagonal structure of the 
circumpolar jet is that the central vortex is close enough to sweep the inner sec- 
ondary vortices appearing due to development of the barotropic instability of the 
jet, and thus damp the meandering observed in the “jet-only” case. 

The above-described results correspond to initialization with a single dominant 
mode 1 = 6. We carried out simulations initialized with a spectrum of random- 
phase modes with different l = 2, 3, ..., 8, which also leads to a hexagonal jet with 
similar properties as in the reference case of single-mode initialization, although 
it takes a longer time to settle. Fig. 16 illustrates the emergence of l = 6 as the 
dominant mode in this case. 
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Figure 16. Evolution of the logarithms of the normalized amplitudes of the Fourier modes 
of the azimuthal velocity in a simulation of “jet+vortex” configuration in the northern 
hemisphere initialized with modes l = 2, 3, ...,8 with the same (small) amplitudes and 
random phases. V (r) is the same as Fig.l. 


5 Summary and discussion 


In the present paper we used the RSW model and the state-of-the-art Cassini ob- 
servations of the polar hexagonal jet and central vortex (Antunano et al., 2015) for 
a search of dynamical mechanisms that could explain the existence and longevity 
of the Saturn’s northern polar hexagon and the absence of its counterpart at the 
South pole. Insufficient knowledge of the vertical structure of the Saturn’s atmo- 
sphere justifies the use of such vertically integrated model, where all information 
about the vertical structure is condensed in a single parameter, the Rossby defor- 
mation radius R 4 . We used a fit of the observed velocity distributions to perform 
detailed linear stability analyses and to initialize nonlinear RSW simulations. In 
contrast with existing literature we analyzed not only a “jet-only”, but also a 
“jet+vortex” configuration, and showed the importance of including the central 
vortex in the system. 

Our main results are as follows 

(1) There indeed exists a region of parameters where the most unstable mode of 
the jet-only configuration has azimuthal wavenumber l = 6, although small 
changes of parameters may shift this value up or down. For the distance 
actually observed between the northern polar jet and the North Pole, the 
azimuthal wavenumber l = 6 is obtained for a Rossby radius of deforma- 
tion Rd « 3200 km. 

(2) In the Southern hemisphere configuration, the most unstable mode has an 
azimuthal wavenumber higher than 8, with a precise value depending on R ( j. 


21 


This is related to the distance of the southern polar jet from the South pole 
being larger than in the Northern hemisphere. 

(3) The preceding results hold for both “jet-only” and “jet+vortex” cases. The 
difference between the two configurations is revealed by high-resolution fully 
nonlinear simulations initialized either with the most unstable mode of small 
amplitude, or with a random combination of unstable modes, superimposed 
onto the background velocity profile. In the “jet-only” configuration, nonlinear 
saturation of the instability leads to the formation of counter-jets and strong 
distortion of the hexagon at the late stages of the evolution, consistently 
with the standard scenario of saturation of the barotropic instability. On the 
contrary, in the “jet+vortex” configuration, the central vortex stabilizes the 
hexagon, leading at the late stages of the evolution to a configuration close 
to the observed one. 

The persistence of the hexagonal structure in our RSW simulations suggests that 
the combination of central vortex and hexagonal jet structure could be an exact 
solution to the system. Proving this hypothesis is a challenge. Such proof would 
close the problem because, however robust and long-lived is this structure in our 
direct numerical simulations, although these latter are very long by the standards 
of computational fluid dynamics (more than a thousand of inertial periods), they 
remain far short of the time of observations of the Saturn’s hexagon. 

Compared to the existing literature, although the model and the background flow 
configurations are different, our “jet-only” configuration yields a “vortex-street” 
unstable northern polar hexagon similar to the results of Morales-Juberias et al. 
(2011), and the “jet+vortex” configuration yields a “meandering jet” stable north- 
ern polar hexagon, similar to Morales-Juberias et al. (2015), featuring the sharp 
straight-line PV gradients on the six sides of the hexagon. This means that the pres- 
ence of the central vortex could stabilize the northern polar hexagonal jet already 
in a simple barotropic model, while a (baroclinic) vertical shear was necessary in 
Morales-Juberias et al. (2015) for stabilization of the hexagonal “meandering jet” . 


In spite of successful representation of the hexagon and explanation of its ab- 
sence at the opposite pole, we still have two caveats in our scenario. The first one 
is too high phase velocity of the hexagon. We present in Table 1 the values of 
the angular (or phase) velocity c of the hexagon for different values of parame- 
ters, as follows from the linear stability analysis of both the “jet-only” and the 
“jet+vortex” configurations (it should be stressed that the angular velocities re- 
trieved from the direct nonlinear RSW simulations are close to results of the linear 
stability analysis). In both cases, the hexagon in our model is not quasi-stationary 
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Table 1 


Phase velocity of the hexagon 


6 

R d [km] 

c [r ad. /o], sole jet 

c [rad./o], jet+vortex 

0.08 

3200 

7.33 x 10 -3 

7.83 x 10 -3 

0.09 

3200 

8.40 x 10" 3 

8.88 x 10" 3 

0.10 

3200 

9.38 x 10 -3 

9.88 x 10“ 3 

0.11 

2400 

7.33 x 10 -3 

7.51 x 10“ 3 

0.12 

2400 

7.71 x 10" 3 

8.20 x 10" 3 

0.14 

2400 

9.06 x 10“ 3 

9.56 x 10“ 3 


as in the observations (Sanchez-Lavega et al., 2014) (a; = +0.0128 ± 0.0013 de- 
gree longitude per Saturn day), but is rotating with c ~ +8 x 10 -3 radian. / 0 , 
i.e. uj = +6 degree longitude per Saturn day, meaning that the hexagonal pattern 
achieves a complete rotation in Saturn’s rotating frame in 60 days, which is much 
higher than in observations. Sanchez-Lavega et al. (2014) use, however, Voyager’s 
System III reference frame based on Saturn’s Kilometric Radiation (SKR). Not 
only Cassini’s SKR measurements yield a 8— min longer rotation period with re- 
spect to Voyager’s (Gurnett et al., 2007), but alternative methods based on gravity 
measurements (Helled et al., 2015), or potential vorticity mapping (Read et al., 
2009b), yield a 5- min to 7-min smaller rotation period, with respect to Voyager’s. 
Assuming our shallow-water simulations are consistent with the atmospheric dy- 
namics framework adopted in Read et al. (2009b), which entails that our simulated 
+6°/day drift rate would translate to a slower +3°/day drift rate in the Voyager 
SKR reference, and even a drift rate close to 0°/day in the Cassini SKR reference. 


The second one is the exaggerated zonal velocity in low latitudes resulting from 
the saturation of the instability. The second caveat is most probably due to the 
limitations of the gamma-plane model, simulations with the RSW model on the 
full sphere would be sufficient to answer the question whether this is a defect of 
the model, or of the gamma-plane approximation. The first one can be probably 
resolved by introducing a vertical structure in the model, such as an active lower 
layer, although this would require a number of ad-hoc hypotheses which we are 
reluctant to introduce. We should recall that baroclinic waves are slower than the 
barotropic ones, so a baroclinic structure would necessary slow down the system, 
the quantitative estimates depending on the details of the vertical structure. 
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It should be also emphasized that we were working with a dissipationless version 
of the model. One of the reasons was to reduce the number of ad hoc hypotheses. 
Of course, an explicit (turbulent) viscosity and diffusion terms could be added in 
the r.h.s. of momentum and mass conservation equations of the model. In order to 
maintain any non-trivial steady state a forcing should be also added in this case 
to prevent dissipative decay. Both forcing and dissipation (the values of turbu- 
lent viscosity and diffusivity) would be necessarily ad hoc and should be adjusted 
to produce the observed mean flow. Above, we were pursuing the idea that the 
Saturn’s polar hexagon results from an instability of this mean flow, i. e. the po- 
lar jet plus polar vortex system. The only difference, if this system is considered 
as a solution of forced-dissipative equations instead of non-forced non-dissipative 
ones, is that eigen-frequencies of its small perturbations would be subject to ad- 
ditional viscous damping, which is necessarily weak in view of longevity of the 
structure in question. A radiative damping can be also introduced in the system. 
In shallow-water modeling of the atmosphere such damping is usually represented 
by relaxation of the thickness of the layer towards some radiative equilibrium 
value with some characteristic relaxation time. As this time is presumably very 
long for Saturn, it was neglected. Another reason of working with inviscid model 
is that numerical scheme used for simulations of the saturation of the instability 
is practically dissipationless, which allows for very long-time runs. 


A Influence of moist convection upon the devloping instability 


Radiative damping mentioned above is an example of diabatic flux in the system. 
Another one, which was recently invoked for explanation of the emergence of Sat- 
urn’s polar vortex with the help of the RSW model (O’Neill et al., 2016), and is 
recurrently evoked to explain the maintenance of alternating jets on giant planets, 
cf. e.g. (Showman, 2007) is moist convection. In relation to the discussion above, 
the moist waves are slower than the “dry” ones, so the moist effects could, in 
principle, explain the observed slowness of the hexagon’s rotation. The condensa- 
tion and moist convection can be included in the RSW model in a self-consistent 
way (Bouchut et al., 2009). For this, a new quantity, the bulk moisture content in 
the atmospheric column Q is introduced. It is conserved in the absence of phase 
transitions, but acquires a sink if condensation is switched on. The condensation 
happens whenever the threshold of saturation Q s is crossed. An extra convective 
flux through the upper surface of the atmospheric layer is added and linked to the 
latent heat of condensation. The resulting equations of so-called moist-convective 
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Figure 1. Snapshots of humidity initialized with uniform through the whole domain value 
qo = 0.89 corresponding to the “jet-h vortex” (left and middle panel at t — 150, lOOO/o -1 ,) 
and u jet-only ” (right panel at t — 150 fo -1 ) configurations. 

rotating shallow water (mcRSW) read: 

/ 

d t v + (v-'V)v + fz x v = —gVh, 

< dth + V- (vh) = -PP, (A.l) 

K d t Q + V-(Qv) = -P, 

with a coefficient f3 related to background stratification, and a simple relaxation 
parameterization for condensation P with relaxation time r: 

P= Q ~ Q * H(Q-Q a ). (A. 2) 

T 

These new elements can be easily incorporated in the finite-volume numerical 
scheme we are using. 

Before making simulations with the full-strength mcRSW, we can get useful in- 
sights from analysing the evolution of Q in the absence of the condensation sink. 
Q then obeys the conservation equation, the third equation in (A.l) with P = 0. 
The results of previous RSW simulations with added passive tracer Q evolving 
from initially uniform distribution are presented in Fig. 1, for both “jet-only” and 
“jet-|-' vortex” configurations. Fig. 1 shows that the excess of humidity resulting 
from the evolution of the system is concentrated in the vicinity of the jet, and 
does not affect the central vortex. It is generally known, and confirmed by the 
simulations with moist-convective RSW in Lambaerts et al. (2011); Rostami and 
Zeitlin (2017), that condensation and related latent heat release lead to intensifi- 
cation of the cyclones and damping of the anticyclones So we can infer that in the 
absence of additional sources of moisture, condensation would not help to main- 
tain the central vortex, and thus to increase the longevity of the system, but would 
only impact the hexagonal jet on its equatorward flank. Simulations with the full 
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Figure 2. Snapshots of condensation (thick black contours) and isobars 0.6, 0.7 and 0.75 
(grey contours) in a most-precipitating environment, as follows from the simulations 
with moist-convective RSW of Bouchut et a 1. (2009). Humidity is distributed uniformly 
in the flow domain, with initial value Q = 0.48 with a saturation threshold at Q s = 0.51. 
Other initial or saturation values do not change the evolution qualitatively. 

moist-convective RSW (A.l) confirm this prediction, as follows from Fig. 2 where 
evolution of condensation in time is displayed. Again, this result may change in 
the two-layer version of the model, but its introduction requires more information 
on the vertical structure of the Saturn’s atmosphere. 
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